###Empirical Study of Lin 2024

rm(list = ls())
ptm <- proc.time()
set.seed(12345);
setwd('/Users/zlin30/Dropbox/Empirics/SIQ')


data=read.csv(file.choose(), header = TRUE, sep = ",") ##Data file contains X
Friend=read.csv(file.choose(), header = TRUE, sep = ",") ##Data file contains friendship

n=dim(data)[1]

# tau=0.1
# tau=0.25
# tau=0.5
# tau=0.75
tau=0.9

library(Matrix)

  NF = rowSums(Friend);
  invNF =1/NF;#1 over number of friend
  invNF[which(!is.finite(invNF))]=0;
  
  
  intelligence=data$intelligence
  white=data$white
  female=data$female
  income=exp(data$logincome)
  age=data$age
  gpa=data$gpa
  
  intel1=as.numeric(intelligence<3)
  intel2=as.numeric(intelligence>2 &intelligence<5)
  intel3=as.numeric(intelligence>4)
  
  
  exercise=data$exercise
  x=cbind(age,female,gpa,intel2,intel3,white,income)
  y=exercise>2
  
  bne=rbinom(n,1,0.5);
  xnet0=as.matrix(Friend)%*%as.matrix(bne)*invNF; ##Construct network regressor
  
  
  source("sa.nps.R")
  err=10;
  err.tol=1E-2;
  bne0=rep(0,n);
  npsms.reps = 0;
  tot.reps = 100;
  quiet=FALSE;
  nc=0;
  while(err>err.tol&npsms.reps<tot.reps){
    npsms.reps = npsms.reps + 1; 
    npsms=sa.nps(y,x,xnet0,tau,lbc=-10,ubc=10);
    betahat0=npsms$coef;
    hopt0=as.vector(npsms$sigma);
    xx1=cbind(1,xnet0,x);
    bne1 = as.integer(as.matrix(xx1)%*%as.vector(betahat0)>0);
    # err = norm(matrix(bne1-bne0,n,1),"f")/sqrt(n);
    err=mean(abs(bne1-bne0))
    xnet0=as.matrix(Friend)%*%as.matrix(bne1)*invNF;
    bne0=bne1;
    
    if(!quiet) {	
      cat("NPS Estimator: ",npsms.reps," iterations completed \n")
      cat(" Current Error: ",err,"\n \n");
    }
    if(npsms.reps==tot.reps & err>err.tol) {
      cat("The NPSMS algorithm did not converge"); 
    }
  }
  cat("Number of iterations to converge:",npsms.reps,"\n\n")
  sms=as.vector(t(npsms$coef))
  norm=sqrt(sum(sms[2:9]^2))
  est=sms/norm
round(est[2:9],digits=3)

##Resampling
scode=unique(data$scid)
ns=length(scode)
neach=as.numeric(table(data$scid))
scid=data$scid

nres=1000
resamest=matrix(rep(0,8*nres),nres,8)
for(ll in 1:nres){
  cat(ll," resampling\n")
resam=sample(scode,ns,replace=TRUE)

subdata=data[scid==resam[1],]
for(i in 2:ns){
  subdata=rbind(subdata,data[scid==resam[i],])
}

subF=Friend[scid==resam[1],scid==resam[1]]
for(i in 2:ns){
  subF=bdiag(as.matrix(subF),as.matrix(Friend[scid==resam[i],scid==resam[i]]))
}

nr=dim(subdata)[1]
NFr = rowSums(subF);
invNFr =1/NFr;#1 over number of friend
invNFr[which(!is.finite(invNFr))]=0;


intelligencer=subdata$intelligence
whiter=subdata$white
femaler=subdata$female
incomer=exp(subdata$logincome)
ager=subdata$age
gpar=subdata$gpa

intel1r=as.numeric(intelligencer<3)
intel2r=as.numeric(intelligencer>2 &intelligencer<5)
intel3r=as.numeric(intelligencer>4)


exerciser=subdata$exercise
xr=cbind(ager,femaler,gpar,intel2r,intel3r,whiter,incomer)
yr=exerciser>2

bner=rbinom(nr,1,0.5);
xnet0r=as.matrix(subF)%*%as.matrix(bner)*invNFr; ##Construct network regressor


source("sa.nps.R")
errr=10;
err.tolr=1E-2;
bne0r=rep(0,nr);
npsms.repsr = 0;
tot.repsr = 100;
quietr=FALSE;
ncr=0;
while(errr>err.tolr&npsms.repsr<tot.repsr){
  npsms.repsr = npsms.repsr + 1; 
  npsmsr=sa.nps(yr,xr,xnet0r,tau,lbc=-10,ubc=10);
  betahat0r=npsmsr$coef;
  # hopt0r=as.vector(npsms$sigma);
  xx1r=cbind(1,xnet0r,xr);
  bne1r = as.integer(as.matrix(xx1r)%*%as.vector(betahat0r)>0);
  # err = norm(matrix(bne1-bne0,n,1),"f")/sqrt(n);
  errr=mean(abs(bne1r-bne0r))
  xnet0r=as.matrix(subF)%*%as.matrix(bne1r)*invNFr;
  bne0r=bne1r;
  
  # if(!quietr) {	
  #   cat("NPS Estimator: ",npsms.repsr," iterations completed \n")
  #   cat(" Current Error: ",errr,"\n \n");
  # }
  # if(npsms.repsr==tot.repsr & errr>err.tolr) {
  #   cat("The NPSMS algorithm did not converge"); 
  # }
}
# cat("Number of iterations to converge:",npsms.repsr,"\n\n")
smsr=as.vector(t(npsmsr$coef))
normr=sqrt(sum(smsr[2:9]^2))
estr=smsr/normr
resamest[ll,]=round(estr[2:9],digits=3)
}


est1=resamest[,1]
est2=resamest[,2]
est3=resamest[,3]
est4=resamest[,4]
est5=resamest[,5]
est6=resamest[,6]
est7=resamest[,7]
est8=resamest[,8]

ci1=quantile(est1,c(0.025,0.975))
ci2=quantile(est2,c(0.025,0.975))
ci3=quantile(est3,c(0.025,0.975))
ci4=quantile(est4,c(0.025,0.975))
ci5=quantile(est5,c(0.025,0.975))
ci6=quantile(est6,c(0.025,0.975))
ci7=quantile(est7,c(0.025,0.975))
ci8=quantile(est8,c(0.025,0.975))



ci=t(cbind(ci1,ci2,ci3,ci4,ci5,ci6,ci7,ci8))
round(cbind(est[2:9],ci),digits=3)
# neachres=rep(0,ns)
# for(i in 1:ns){
# neachres[i]=neach[scode==resam[i]]
# }
# 
# sum(neachres)==dim(subdata)[1]



proc.time() - ptm
